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Abstract 

The Lyapunov spectrum describes the exponential growth, or decay, of infinites- 
imal phase-space perturbations. The perturbation associated with the maximum 
Lyapunov exponent is strongly localized in space, and only a small fraction of all 
particles contributes to the perturbation growth at any instant of time. This fraction 
converges to zero in the thermodynamic large-particle-number limit. For hard-disk 
and hard-sphere systems the perturbations belonging to the small positive and large 
negative exponents are coherently spread out and form orthogonal periodic struc- 
tures in space, the "Lyapunov modes" . There are two types of mode polarizations, 
transverse and longitudinal. The transverse modes do not propagate, but the longi- 
tudinal modes do with a speed about one third of the sound speed. We characterize 
the symmetry and the degeneracy of the modes. In the thermodynamic limit the 
Lyapunov spectrum has a diverging slope near the intersection with the abscissa. 
No positive lower bound exists for the positive exponents. The mode amplitude 
scales with the inverse square root of the particle number as expected from the 
normalization of the perturbation vectors. 
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1 Introduction 



Recently, the concepts of dynamical systems theory and molecular dynamics 
simulations have been applied to the phase-space dynamics of models repre- 
senting simple fluids and solids [1,2]. Due to the convex and dispersive na- 
ture of the atomic surfaces, the phase trajectory of such a system is highly 
Lyapunov unstable resulting in exponential growth, or decay, of small (in- 
finitesimal) initial perturbations along specified directions in phase space. The 
respective rate constants, the Lyapunov exponents, are taken to be ordered 
by size with Ai > representing the maximum exponent. The whole set, 
{A;, / = 1, . . . , 2dN}, is known as the Lyapunov spectrum, where d denotes the 
dimension of space, and N is the number of particles. For fluids in nonequi- 
librium steady states close links exist between the Lyapunov spectrum and 
various macroscopic dynamical properties - transport coefficients, irreversible 
entropy production, and the Second Law of thermodynamics [1,2,3,4,5,6]. This 
important result provides the motivation for us to examine also the spatial 
structure of the perturbations associated with the various exponents. Here we 
restrict ourselves to systems in thermodynamic equilibrium. Earlier accounts 
of our work have been published in Refs. [7,8,9,10,11,12,13,14]. 

Since the pioneering work of Bernal [15] with steel balls, and the seminal 
computer simulations of Alder and coworkers [16,17], hard ball systems are 
considered good models for the structure of "real" dense fluids. They serve as 
reference systems for highly-successful perturbation theories of liquids [18,19]. 
The dynamics, however, is rather different, and it is rewarding to study this 
difference, not only with the established machinery of correlation functions, 
but also from the viewpoint of dynamical systems theory. In this paper we 
discuss our results on the Lyapunov spectra and the space dependence of the 
associated perturbation vectors for two-dimensional hard disk systems in equi- 
librium. For comparison we also present related results for two-dimensional 
soft-disk fluids. 

This work has three main objectives: 

1) In Section 3 we demonstrate that the perturbation associated with the 
maximum Lyapunov exponent, Ai, is strongly localized in space: only a small 
fraction of all particles contributes actively to the growth of the perturbation 
norm at any instant of time. This localization persists in the thermodynamic 
limit, and, even more surprisingly, the fraction of simultaneously-contributing 
degrees of freedom converges to zero. 

2) For small exponents close to the intersection with the abscissa the Lya- 
punov spectrum has a step-like appearance as is depicted in Fig. 1 for a 
two-dimensional 1024-disk fluid with particle density p = 0.7. The steps are 
caused by degenerate exponents, and the associated perturbations form coher- 
ent wave-like patterns in space with well-defined wave vectors and polariza- 
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Fig. 1. Lyapunov spectrum of a 1024-disk system in a square simulation box with 
periodic boundaries. The density p = 0.7. A reduced index 1/2N is used on the 
abscissa. The steps are magnified in the inset. The non-normalized index, I, is used 
there. The spectrum is defined only for integer values of I. 

tion. We refer to these patterns as "Lyapunov modes" . We have found them in 
hard-disk and hard-sphere systems in two and three dimensions [7,8,9,10], re- 
spectively, and in hard-dumbbell systems modelling fluids consisting of linear 
molecules [11,12,13,14]. In Section 4 we characterize the symmetry and polar- 
ization of these modes, where the emphasis is again on the thermodynamic 
limit for bulk systems: N and V — > oo such that p =const. V is the area of the 
box. We show that in this limit the Lyapunov spectrum has a vertical slope 
at the intersection point with the abscissa. No gap, or lower positive bound, 
exists for the positive exponents, and the smallest positive exponent converges 
linearly to zero with inverse system size. 

3) How do the modes depend on the intermolecular potential? In Section 6 
we discuss simulations of two-dimensional particle systems interacting with a 
continuous repulsive Weeks-Chandler- Anderson (WCA) potential. They reveal 
that the modes are extremely unstable and not easily observed. We relate this 
phenomenon to the strong fluctuations of the time-dependent local Lyapunov 
exponents. 

Since our discovery of the Lyapunov modes the theoretical foundation of the 
Lyapunov modes has been examined in some detail [20,21,22,23,24,25,26]. This 
work will be summarized in the concluding Section 7. 
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To facilitate the interpretation of our numerical results we outline our numer- 
ical methods for the computation of Lyapunov spectra for hard-disk systems 
[27,7] in the following section. 



2 Lyapunov spectra and numerical procedure 

The instantaneous state of a planar particle system is given by the AN- 
dimensional phase-space vector T = {q iy p^, ; i — 1, . . . , N}, where = (x iy i/i) 
and pj = (p x ,i,P y ,i) denote the position and linear momentum of particle i, 
respectively. An infinitesimal perturbation or offset vector 5F = {5q«, Spf, i = 
1, . . . , N} evolves according to motion equations obtained by linearizing the 
dynamical equations for T(t). For hard disks this amounts to a linearization of 
the collision map, which maps pre-collision into post-collision states in phase 
space [27]. There exist AN orthonormal initial vectors (<5r z (0); 1 = 1,..., AN} 
in tangent space, such that the Lyapunov exponents 

exist and are independent of the initial state [29,30]. Geometrically, the Lya- 
punov spectrum describes the stretching and contraction along linearly-in- 
dependent phase space directions of an infinitesimal hypersphere co-moving 
with the flow. According to the conjugate pairing rule [31,32,33] for sym- 
plectic systems, the exponents appear in pairs, and the pair sum vanishes, 
A; + A 4 at + i_z = 0. Thus, only the positive half of the spectrum {\\<i<2n} 
needs to be calculated. Furthermore, six of the exponents, {A2at-2</<2JV+3}, 
vanish as a consequence of the conservation of energy, momentum, and center 
of mass, and of the non-exponential time evolution of a perturbation vector 
parallel to the phase flow. 

For our numerical work we use a variant [27,7] of the classical algorithm by 
Benettin et al. and Shimada and Nagashima [34,35], where the reference tra- 
jectory and the orthonormal set of perturbation vectors are simultaneously 
evolved in time. The latter are periodically re-orthonormalized with a Gram- 
Schmidt (GS) procedure after consecutive time intervals At GS , which were 
chosen according to Atcs — r 2/5. Here, T2 = {(3m) 1 ^ 2 / (27i 1 ^ 2 pa) is the colli- 
sion time per particle, (3 = (/cbT) -1 , T is the temperature, k B is Boltzmann's 
constant, and a is the disk diameter. The Lyapunov exponents are determined 
from the time-averaged renormalization factors. 

Fig. 1 shows the positive branch of the Lyapunov spectrum for a 1024-disk 
system in a square simulation box with periodic boundaries. If not stated 
otherwise, reduced units are used for which the particle diameter a, the particle 
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Fig. 2. Squared norm of the individual particle perturbations, |<5qj| 2 + |<5pj| 2 , for 
the offset vector associated with Ai, plotted at the positions q« = (xi,yi) of the 
particles. The ensuing surface is linearly interpolated over a periodic grid covering 
the simulation box. The system consists of 1024 hard disks in a square box with 
periodic boundaries, at a density p = 0.7, and a temperature T = K/N = 1. 

mass m, and the kinetic energy per particle, K/N = ksT, are unity. Thus, the 
unit of time is t* = [ma 2 N/ K) 1 / 2 , and the Lyapunov exponents are given in 
units of 1/t*. The simulation box has lengths L x and L y . In all our examples 
the density p = N/(L x L y ) is 0.7, which is slightly smaller than the fluid-to- 
solid phase transition density. 



3 Localized mode for the maximum exponent 

The maximum exponent is dominated by the fastest dynamical events. There is 
strong numerical evidence that the thermodynamic limit {N, V — > oo, N/Vcon- 
stant} for Ai and, hence, for the whole spectrum exists [27,13]. The associated 
perturbation is strongly localized in space [28,12,13]. This is demonstrated by 
projecting the perturbation vector 5Yi onto the four-dimensional subspaces 
spanned by the components belonging to the individual particles. The squared 
norm of this projection, \5(\/\ 2 + |5pi| 2 , indicates how active a particle % is en- 
gaged in the growth process characterized by Ai [13,12]. In Fig. 2 it is plotted 
for a 1024-disk system at the particle positions (x;,?/;), where the ensuing 
surface is interpolated over a periodic grid covering the simulation box. The 
figure refers to an instantaneous configuration of a well-relaxed system. The 
fastest expansion activity is clearly confined to a few very active zones. 

To understand this localization we note that the offset-vector dynamics is 
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Fig. 3. Particle-number dependence of the localization measure for the perturba- 
tion belonging to the maximum exponent, Ci 5 e- The threshold is 0.98. The line 
represents a power law, 25iV~ ' 76 . 

governed by linear equations and that a perturbation component has only a 
fair chance of further growth if its value was already above average before a 
collision. Each global renormalization step tends to reduce the (already small) 
components of the other non-colliding particles even further. There is a com- 
petition for growth which favors the particles with the largest perturbation 
components. The few active zones move diffusively in space and perform oc- 
casional jumps. 

The localization persists in the thermodynamic limit. To show this [13] we 
order all squared components [51^; j = 1, . . . , 4N of a perturbation vector STi 
according to size. By adding them up, starting with the largest, we determine 
the smallest number of terms, n, which are required for the sum to exceed a 
threshold 9. Q ; e = n/AN is a relative measure for the number of components 
contributing to A;: 



Obviously, Q ; i = 1. In Fig. 3 the localization measure, for / = 1 and for 
a threshold 6 = 0.98, is shown as a function of the particle number N. It 
obeys a power law with a negative power, and converges to zero for iV — > oo. 
The ratio of tangent- vector components (and, hence, of particles) contributing 
significantly to the maximum Lyapunov exponent vanishes in that limit. A 
power-law dependence for the localization measure is a very general result, 
although the parameters are sensitive to the interaction potential [9,36]. 
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Fig. 4. Dependence of the localization measure C;,e=o.98 on the reduced Lyapunov 
index 1/2N. The system consists of N = 1024 particles, and the density p is 0.7. 
The inset provides an enlargement of the mode regime. Here, the Lyapunov index 
I is not normalized. 

If the localization measure is plotted for all positive exponents with index 
1 < I < 2N as depicted in Fig. 4, one observes a fast increase for small /, 
where the full dot indicates the result for / = 1. This is an indication that 
the localization does not persist for large / and that the perturbation growth 
starts to involve more and more particle components. The inset shows the 
interesting regime for which modes have been observed. 



4 Perturbation vectors for the six vanishing exponents 

To be explicit, we arrange the components of the 4iV-dimensional state vector 
according to 

r = . . .,x N ,y N ;p Xt i,p Vt i, . . . ,p x ,n,P v ,n) ■ (3) 

During the streaming motion in between instantaneous elastic collisions no 
forces act between the particles, and a normalized vector parallel to the phase 
flow r becomes 
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(4) 



A perturbation vector parallel to e x does not grow exponentially with time 
and contributes one vanishing exponent to the full spectrum [5]. 

The motion in the extended 4iV-dimensional phase space is constrained to hy- 
persurfaces generated by center-of-mass conservation, (J2 x f=i — const, Y^VjLi 
= const), by momentum conservation, (J2jLiPx,j = 0, Y^jLiPyj = 0), and by 
the conservation of energy (Y^f=i(Px,j ~^Py,j)/^ = ^0- -^U respective vectors 
orthogonal to these surfaces also denote directions of non-exponential phase- 
space expansion or compression: 



e 2 = (l/ v / iV)(l,0,...,l,0;0,0,...,0,0) (5) 

e 3 = (l/v / iV)(0,l,...,0,l;0,0,...,0,0) (6) 

e 4 = (l/ v / iV)(0,0,...,0,0;l,0,...,l,0) (7) 

e 5 = (l/v / iV)(0,0,...,0,0;0,l,...,0,l) (8) 

e 6 = (1/ V2K) (0,0,..., 0, 0; p Xtl ,p Vtl , Px,n,Pv,n) • (9) 



The vectors ei to e$ are orthonormal and form the basis of an invariant sub- 
space, the central manifold, which contains also the six perturbation vectors 
{ST i, 2N — 2 < I < 2n + 3} responsible for the vanishing exponents. An 
expansion in the natural basis gives 

6 

<$r, = , / = 2iV-2,...,2iV + 3 . (10) 

i=\ 

The matrix of the projection cosines, ol\^ = (STi ■ e,), obeys the normalization 
conditions 

6 2N+3 

EW 2 = i , E («w) 2 = i, 

i=l l=2N-2 

if the total momentum vanishes and if K = N as is always the case in our simu- 
lations. As an example we list such a matrix in Table 1 for a 400-disk fluid. We 
note that these numbers are the asymptotic values for a well-relaxed system 
and do not change with time. They always appear in block-diagonal form, but 
other than that they depend on the initial conditions for the simulation. The 
block-diagonal structure indicates that the subspaces C + =span{ei, e 2 , e 3 } and 
C~ =span{e 4 , e 5 , e 6 } are asymptotically ordered by the Gram-Schmidt proce- 
dure due to linear perturbation growth/decay during the streaming motion 
between successive collisions. This may be seen explicitly from the linear mo- 
tion equations for intercollisional streaming and the linearized collision map. 
For example, perturbations ST, parallel to ei, e 2 , or e 3 at t — immediately 
after a collision, will continue to be so at a later time t > 0. However, pertur- 
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Table 1 

Projection cosines of the perturbation vectors in the central manifold onto the nat- 
ural basis vectors for a well-relaxed configuration of a 400-disk fluid with density 
p = 0.7 and a total kinetic energy K = N . 
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(*l,i 


ai,2 


«Z,3 




«;,5 


"Z,6 


398 


0.431 


0.297 


-0.852 











399 


0.549 


-0.836 


-0.013 











400 


-0.716 


-0.462 


-0.523 











401 











-0.138 


-0.607 


-0.783 


402 











-0.620 


0.669 


-0.409 


403 











-0.772 


-0.429 


0.469 



bations ST parallel to e 4 , e 5 , and e 6 at time zero, evolve according to e 2 t + e 4 , 
e 3 t + e 5 , and e±t + e 6 , respectively, for t > previous to the next collision. 



5 Lyapunov modes 

5. 1 Classification of the modes 

Next, we study the non- localized perturbations associated with the small posi- 
tive exponents. The simulation box is a square with L x = L y = L. The inset in 
Fig. 1 gives an enlargement of the step-like structure for a 1024-disk spectrum. 

The first step consists of four degenerate exponents, for which I = 2045, 2044, 
2043 and 2042. On the left-hand side of Fig. 5 we plot the positional pertur- 
bation components, Sxi and 

5yi, of the particles at their positions (xi,yi) in the simulation box, where 
the surfaces have been linearly interpolated on a periodic grid covering the 
simulation box. One finds periodic patterns like sines and cosines, with a phase 
shift of it/2 between them. The wave vectors for Sx and Sy are orthogonal to 
the y and x axes, respectively, and we refer to these modes as transverse. Their 
wave number, k = |k|, is the smallest consistent with the periodic boundaries, 
ki = 2n/L for 2042 < I < 2044. The same transverse modes, in another 
representation, are shown on the right-hand side of Fig. 5, where the particle 
perturbations (Sxi, Syi) for all particles i are plotted as small arrows at the 
respective particle positions. Transverse modes do not propagate. 

The next-larger exponents responsible for the second step in the inset of Fig. 1 
have multiplicity eight, 2034 < I < 2041. Their modes are called longitudinal, 
since the perturbations 5x and Sy are parallel to their respective wave vectors 
with wave number 2n/L. Longitudinal modes are propagating, with a phase 
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Fig. 5. Transverse modes with the smallest-possible wave number for I = 2045, 2044, 
2043 and 2042. The fluid consists of 1024 disks with a density p = 0.7. See the main 
text for details. 
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Fig. 6. Dependence of the small positive Lyapunov exponents on the wave number 
k for various systems containing between 100 and 1024 particles. L refers to longi- 
tudinal, and T to transverse modes. The density of the hard-disk fluid is 0.7, the 
kinetic energy per particle, K/N, is unity. 

velocity roughly equal to one third of the sound speed [8] . 

Continuing up the steps, one encounters transverse modes for 2030 < / < 
2033 with multiplicity 4 and wave vectors pointing along the diagonals of 
the simulation box, k = 2y / 27r/L, followed by longitudinal modes for 2030 < 
/ < 2033, with multiplicity 8 and wave vectors also along the diagonals, k = 
2 v / 27r/L, and so on. In the general case of a rectangular simulation box, all 
modes look like waves with n x nodes in x direction and n y nodes in y direction, 
such that 



k = 2tt 



\ 



+ 



n,. 



'ID 



The case n x = n y = corresponds to perturbation vectors in the central ma- 
nifold and has been treated in the previous section. For all positive exponents 
the momentum perturbations Sp x and Sp y are strictly parallel to the respective 
positional perturbations 5x and 5y, and are therefore not included in Fig. 5. 

All modes may be classified either as transverse or longitudinal, where the 
dependence of the Lyapunov exponents on k differs for each class. This is 
shown in Fig. 6, in which all small exponents, plotted as a function of k, lie on 
two curves, one for the transverse (T) and one for the longitudinal (L) modes. 
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Fig. 7. Reciprocal lattice with point spacing equal to 2-ir/L, for sine-like patterns 
(open circles) and cosine- like patterns (dots). Transverse modes do not require a 
sense of direction, and only one half plane is needed (left), whereas for longitudinal 
modes the full reciprocal lattice is required (right). The circle connects allowed wave 
vectors with k = 2y / 5vr/L. 



To second order in k, these "dispersion relations" may be approximated by 



A T = 2.48(3)A; + 0.23(7)A; 2 + C(A; 3 ) (12) 
A L = 3.13(7)A; + 1.2{2)k 2 + 0(k 3 ) . (13) 

The standard deviations affect the last digit of the fit parameters and are given 
in brackets. These second-order approximations are indicated as dashed lines 
in Fig. 6. It has been argued by Mareschal and McNamara [23] that the k 2 
term is negligible. However, the data, in particular the longitudinal exponents, 
are better represented if this term is kept, although it is unessential for the 
discussion below. 



5.2 Mode degeneracy 



The degeneracy of the Lyapunov exponents is a consequence of the square 
symmetry of the periodic simulation box and the phase shift of 7r/2 between 
sine and cosine-like patterns, allowing for different orthogonal perturbations 
with the same k. For the non-propagating transverse modes the sense of di- 
rection is not needed. It suffices to consider all wave vectors in a half plane 
of the reciprocal lattice, as on the left-hand side of Fig. 7. The full and 
open circles refer to sine and cosine-like waves, respectively, which have to be 
counted separately. Degenerate modes belong to all lattice points on a circle 
with a specified radius k. Their multiplicity is denoted by M k . For the propa- 
gating longitudinal modes the full reciprocal lattice plane is needed as on the 
right-hand side of Fig. 7. As a consequence, the degeneracy of the longitudinal 
modes is always twice that of the transverse modes for a given k. With the 
circles in Fig. 7 we select degenerate modes with k = 2^/5ir/L as examples. 
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Fig. 8. Small positive Lyapunov exponents for a 1024-disk system: the crosses denote 
simulation results, the smooth curve indicates the reconstruction to second order in 
k as described in the main text. 

The transverse modes have multiplicity 8 and appear for I = 2010, . . . , 2017 
in the 1024-disk spectrum of Fig 8. The longitudinal modes have mutiplicity 
16 and are indexed by I = 1978, . . . , 1993. Of course, the relative position of a 
degenerate group of exponents in the spectrum is determined by the dispersion 
relations given above. 



5. 3 Mode amplitudes and stability of transverse modes 

Each mode with index / is viewed here as a four-dimensional field on the 
periodic domain of the simulation box, and will be denoted byQ, 

V t (x, y) = {5x(x, y), 5y(x, y), 5p x (x, y), 5p y (x, y)}i . (14) 



Its norm, 

\Vi(x, y )\ = i(5xf + (Sy) 2 + (5p x ) 2 + (5 Py ) 2 ]] /2 , (15) 



Sometimes it is advantageous to label the fields, Vi, not by / but by the wave- 
vector norm k and an integer m distinguishing different modes within a degenerate 
group, 1 < m < Mfc 
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is also a periodic function on the domain with amplitude Pq. Since STi is 
normalized, it follows that 

L L 

J J \Vi(x,y)\ 2 dxdy = l. (16) 
o o 



As an example we consider again the transverse modes for / = 2N— 3 belonging 
to the first step in the spectrum. For iV = 1024 such a mode is shown at the 
top of Fig. 5. For any N it has the simple structure 



V 2 n-3(x, y) = {(5x) sin(% + <f> y ), (5y) sm(kx + <f> x ), 

(5p x ) sin(% + (f> y ), (5p y )o sin(kx + x )}2jv-3 , (17) 

where k = 2n/L. The phases <p x and <p y are equal for this particular example, 
but in a more general case they may differ by n. 

The mode amplitude may be expressed in terms of the component amplitudes 
according to 

Po N ~ 3 = [(Sx)l + (5y)l + (5p x )l + {5 Py )l} 1 / 2 . (18) 



The box-size and, hence, the particle-number dependence of this quantity 
follows immediately from the normalization condition Eq. (16): 

iT' 3 = (^) • (19) 



The amplitude vanishes in the thermodynamic limit [20] . In Fig. 9 the depen- 
dence of Pq N ~ 3 on \fN = L is shown. The dots are numerical results, where 
the mode amplitude is obtained from fits of the components of Eq. (17) to 
the data. The smooth line is the result predicted by Eq. (19). Similar results 
are obtained for the other modes. The good agreement between theory and 
experiment is an indication for the remarkable stability of the modes. 

It should be pointed out that only the mode amplitude derived from the vector 
norm, Eq. (15), displays the correct system-size dependence. The amplitudes 
of the vector components do not. Even (Sx)q and (5y)o do not agree in general. 
The amplitudes of the field components depend on the initial conditions of the 
simulation, but Pq does not. 
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Fig. 9. The dots indicate the 1/VN dependence of the transverse-mode amplitude 
Pq N ~ 3 corresponding to the first nonvanishing exponent A2jv-3. The smooth line is 
the prediction of Eq. (19). 



5.4 Spectrum reconstruction near the thermodynamic limit 



If the dispersion relations are known for a fluid of a given density, the rules 
outlined above allow a reconstruction of the most-interesting spectral range for 
small positive exponents. If the second-order approximations of Eqs. (12) and 
(13) are used, the reconstructed spectrum for the 1024-disk system is indicated 
in Fig. 8 by a smooth line. As expected, it agrees very well with the direct 
simulation results for / > 1962 corresponding to small k. Since Eqs. (12) and 
(13) apply to the k — » limit, very large systems may be studied which are 
not accessible by direct simulation, now and in the foreseeable future. In Fig. 
10 a reconstructed spectrum for a fluid containing 40, 000 disks with a density 
0.7 is shown, which may be considered to be close to the thermodynamic limit. 
Only a very small part of the spectrum is actually shown, for which k does 
not exceed 0.3. We may fit a power law to these spectral points, 



A; = a 



1 2N 



' °-" < 2iV < 1 ' (20) 
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Fig. 10. Part of a reconstructed Lyapunov spectrum for a fluid of density 0.7, con- 
sisting of 40,000 disks, is shown by the dashed line. The smooth line is a fit of these 
points to a power law, Eq. (20). For comparison, the spectrum for 1024 disks is also 
shown by the points. On the abscissa a normalized index 1/2N is used. 

with a = 7.81 ± 0.04 and (3 = 0.52 ± O.O1.0 The slope of the spectrum near 
the intersection with the abscissa, 



d\. 



lim — — r ~ lim 

l/2N^l d(l/2N) l/1N-*l 



I 

2N 



-0.48 



(21) 



diverges for 1/2N — > 1. Since the small positive exponents approach zero 
linearly with k as expressed by the dispersion relations above, no gap appears 
in the spectrum. No lower positive bound exists for the positive Lyapunov 
exponents. An analogous statement for the negative exponents of the spectrum 
holds due to the conjugate pairing symmetry of the spectrum. 

We have also extended these considerations to quasi-one-dimensional systems 
[36,37], for which the aspect ratio of the simulation box, A = L y /L x , converges 
to zero in the large-particle limit. This is achieved by fixing the box size L y 
in y direction, typically a few particle diameters, and letting L x become large 
such that the particle density remains constant. The Lyapunov modes are more 
obvious in this case, since the reciprocal lattice consists of equidistant points on 

2 We note that for a simple model of a two-dimensional solid the number of vibra- 
tional modes dl between wave numbers k and k + dk is proportional to k. Integrating 
this relation and assuming linear dispersion relations for the Lyapunov exponents 
we obtain A ~ I 1 / 2 [1]. This power closely agrees with our value for (5. 



16 



the abscissa, allowing only multiplicity two for the transverse modes, and four 
for the longitudinal. The limiting spectral slope near the intersection with the 
abscissa remains finite [36] . Very recently, Taniguchi and Morriss [38] have also 
studied other than periodic boundary conditions for very narrow quasi-one- 
dimensional systems, for which the particles cannot cross. Equivalent modes 
are observed. 

So far we have concentrated on the positive branch of the Lyapunov spectrum, 
for which the components of the momentum perturbations are parallel to the 
respective perturbation components for the positions. This is a consequence 
of the linearity of the motion equations in tangent space and the positivity 
of the exponents. For the negative branch of the spectrum the momentum 
components are antiparallel to the position components of the perturbations. 
All our statements concerning the multiplicity and polarization of the modes 
carry over also to this case. 



6 A puzzle provided by the soft potential 

What happens if the interparticle potential is soft? In Refs. [10,9] a com- 
parison was made between soft and hard-disk systems and, surprisingly, no 
modes could be clearly identified. Here we extend this work by considering 
two-dimensional soft fluids with a Weeks- Chandler- Anderson (WCA) interac- 
tion potential, 



We use units for which the particle mass m, the particle diameter a, and the 
energy parameter e are unity. Furthermore, we consider a thermodynamic state 
with a total energy per particle, E/N, also equal to unity. Since there is a po- 
tential energy in this case, a comparison between soft and hard disks is done for 
equal temperatures T = K/N (and not for equal energy), where Boltzmann's 
constant is also set to unity. All hard-disk Lyapunov exponents quoted in this 



section are therefore rescaled by a factor y (K/N)wca/ \K /N) hd, where K is 
the total kinetic energy. Since the soft-potential spectra are computationally 
more expensive, we restrict our examples to N < 100. 

The spectra for the soft- and hard-disk fluids in Fig. 11 are surprisingly dif- 
ferent, both in shape and in absolute value. Most conspicuously, no steps are 
observed in the WCA case, whereas for hard disks at least the lowest step is 
well developed. This comparison provides us with the following facts [9,36]: 



<j>(r) = < 




(22) 



r > 2 1 / 6 o- . 
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Fig. 11. Lyapunov spectra for 100 soft (WCA) and hard (HD) disks at a density 
p = 0.7 and a temperature T = 0.75. The simulation box is a square with periodic 
boundary conditions 



1) The overall shape of the WCA spectrum may be approximated by a power 
law similar to Eq. (20), but with a power of the order of 0.4, (which depends 
slightly on the fitting range). A similar behavior has been found for three- 
dimensional WCA fluids, and has been interpreted there in terms of a simple 
Debye model for the distribution of vibrational frequencies in solids, anticipat- 
ing the Lyapunov modes. For details we refer to Ref. [1], although the systems 
there are too small for a meaningful comparison with the hard-disk fluids. 

2) If the WCA potential is made progressively steeper, its spectrum converges 
towards that of hard disks, and the step structure starts to re-appear [36]. 
This convergence, however, is slow, in particular for the largest exponent Ai. 

3) The perturbation vector associated with Ai is also localized for soft disks, 
and this localization persists in the thermodynamic limit. An instantaneous 
snapshot of such a perturbation for 100 WCA particles looks very similar to 
Fig. 2. Also the localization measure C^q introduced in Section 3 is found to 
obey a power law with a negative power of N [37]: Ci ; e=o.92 = 0.31Af~ a10 for a 
WCA fluid with p = 0.7 and T = 0.75. Only a few particles, localized in space, 
provide the main contributions to the perturbation vector ST\. However, there 
is a very sensitive dependence of the details on the interparticle potential. 

4) As we suspected from the absence of a step structure, no Lyapunov modes 
could be reliably identified, either visually or by Fourier-transformation tech- 
niques. Although it is possible to see patterns resembling transverse modes for 
the smallest exponent, these patterns are transient and too unstable to survive 
time averaging [10]. This observation is consistent with large fluctuations of 
the local time-dependent exponents \{, from which the Lyapunov exponents 
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Fig. 12. Second moments of the time-dependent local exponents for a planar 16-par- 
ticle system. The labels WCA and HD refer to the Weeks-Chandler-Anderson and 
hard disk models, respectively. The re-orthonormalization time Atcs = 0.075, the 
temperature T = 0.75, and the density p = 0.7. 

are obtained by time averaging: 

1 } 

A, = Jim - / \'i{T(t))dt . (23) 
o 



A';(r) depends on the state T{t) the system occupies in phase space at time 
t and, of course, on a long-enough part of the trajectory terminating at this 
point. They may be estimated from 



where t and t + Atcs refer to times immediately after consecutive Gram- 
Schmidt re-orthonormalization steps. In Fig. 12 the second moments of these 
quantities, {(A'f)}, 1 = 1,... , 4AT}, time averaged along the phase trajectory, 
are shown for 16-particle WCA and hard-disk systems, where a re-orthonor- 
malization interval 0.075 is usecQ I = 1 refers to the maximum, and I = 64 to 
the most-negative exponent. The points for 30 < / < 35 correspond to the 6 
vanishing exponents and are not shown. We infer from this figure that for the 

3 The computation of the second moments (A'f ) for hard disks requires some care. 
For small re-orthonormalization intervals Ai^s, the variance (A' 2 ) — (A'^) 2 varies 
with 1/Atcs and diverges for Atcs — ► 0, but the shape of the fluctuation spectrum 
is hardly affected by the size of Atcs- 
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soft WCA particles the fluctuations of the local exponents become very large 
for the Lyapunov exponents describing relatively-weak instabilities with near- 
zero growth rates, / — > 2N. For the hard disk system, however, the relative 
significance of the fluctuations becomes minimal in this regime. 



7 Outlook 



Soon after the discovery of the Lyapunov modes, Eckmann and Gat [20] were 
the first to provide theoretical arguments for the existence of hydrodynamic- 
like Lyapunov modes. They were derived from the translation invariance of 
the systems, and were based on an evolution matrix in tangent space mod- 
eled as a product of independent random matrices. No accompanying real 
space dynamics exists for this model. Most recently, this theory has been 
improved and made more realistic by Eckmann and Zabey [21]. In another 
approach, McNamara and Mareschal [22] isolate the six hydrodynamic fields 
related to the invariants of the binary particle collisions and the vanishing 
exponents, and use a generalized Enskog theory to derive hydrodynamic evo- 
lution equations for these fields. Their solutions are the Lyapunov modes. In 
a more detailed extension of this work a generalized Boltzmann equation is 
used for the derivation of the hydrodynamic evolution equations [23], which 
restricts the quantitative aspects of the theory to small densities. A related 
approach was taken by de Wijn and van Beijeren [24], who have pointed out 
that the modes should be interpreted as Goldstone modes: The breaking of 
a continuous symmetry - translational invariance in our case - implies the 
existence of hydrodynamic-like "soft" modes and long-ranged correlations. Fi- 
nally, Taniguchi, Dettmann, and Morriss have approached the problem using 
periodic orbit theory [25] and master equations [26]. 

We have mentioned in Section 6 that for small positive exponents the Lya- 
punov spectrum is well described by a power-law close to the thermodynamic 
limit, and that no positive lower bound exists in this case. This has raised 
questions about the numerical accuracy and convergence of our results. With 
respect to the former, we have verified that quadruple precision (128 bits) 
arithmetic gives the same results as the usual double-precision arithmetic. 
To answer the question of convergence, we have perturbed an already well- 
oriented set of orthonormal tangent vectors for a well-relaxed system and 
have ascertained that the unperturbed and perturbed sets of vectors become 
identical again after a certain transient time. The relaxation time the tangent 
vectors require to reach their proper orientations has been studied in Ref. [39]. 

The spatial localization of the fastest-growing perturbation described by the 
maximum exponent is a very general result. It is also found for particle sys- 
tems interacting via soft potentials [9], although the details of the power- law 
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dependence of the localization measure C^e depend sensitively on the inter- 
action potential. Also continuous systems in nonequilibrium stationary states 
exhibiting Rayleigh-Benard convection have this localization property [40]. 
There it is found that the localized active zones contributing most to the per- 
turbation vector for large positive Lyapunov exponents are involved in the 
creation and annihilation of defects of the roll pattern, resulting in the break- 
ing or reconnection of Rayleigh-Benard rolls [40] . 
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